Phase-field-crystal study of grain boundary premelting and shearing in bcc iron 



Ari Adland and Alain Karma 
Department of Physics and Center for Interdisciplinary Research on Complex Systems, 
Northeastern University, Boston, MA 02115 USA 

Robert Spatschek 

Max- Planck- Institut fur Eisenforschung GmbH, D-40237 Diisseldorf, Germany 

Dorel Buta and Mark Asta 
Department of Materials Science and Engineering, 
University of California, Berkeley, California, 94720 USA 
(Dated: November 6, 2012) 

We use the phase-field-crystal (PFC) method to investigate the equilibrium premelting and 
nonequilibrium shearing behaviors of [001] symmetric tilt grain boundaries (GBs) at high homolo- 
gous temperature over the complete range of misorientation < 9 < 90° in classical models of bcc 
Fe. We characterize the dependence of the premelted layer width W as a function of temperature 
and misorientation and compute the thermodynamic disjoining potential whose derivative with re- 
spect to W represents the structural force between crystal-melt interfaces due to the spatial overlap 
of density waves. The disjoining potential is also computed by molecular dynamics (MD) simula- 
tions, for quantitative comparison with PFC simulations, and coarse-grained amplitude equations 
(AE) derived from PFC that provide additional analytical insights. We find that, for GBs over an 
intermediate range of misorientation (# m in < 9 < # max ), W diverges as the melting temperature 
is approached from below, corresponding to a purely repulsive disjoining potential, while for GBs 
outside this range [6 < (9 m i n or # max < 9 < 90°), W remains finite at the melting point, with its 
value corresponding to a shallow attractive minimum of the disjoining potential. The misorientation 
range where W diverges predicted by PFC simulations is much smaller than the range predicted 
by MD simulations when the small dimensionless parameter e of the PFC model, where e -1 / 2 is 
proportional to the ratio of the solid-liquid interface width to the lattice spacing, is matched to 
liquid structure factor properties. However it agrees well with MD simulations with a lower e value 
chosen to match the ratio of bulk modulus and solid-liquid interfacial free-energy, consistent with 
the amplitude-equation prediction that # m i n and 90° — # ma x scale as ~ e 1 ^ 2 . The incorporation of 
thermal fluctuations in PFC is found to have a negligible effect on this range. In response to a shear 
stress parallel to the GB plane, GBs in PFC simulations exhibit coupled motion normal to this 
plane, with a discontinuous change of the coupling factor as a function of 9 that reflects a transition 
between two coupling modes, and/or sliding (shearing of the two grains). Partial sliding is only 
observed over a range of misorientation that is a strongly increasing function of temperature for 
T/Tm > 0.8 and matches roughly the range where W diverges at the melting point. The coupling 
factor for the two coupling modes is in excellent quantitative agreement with previous theoretical 
predictions [J. W. Cahn, Y. Mishin, and A. Suzuki, Acta Mater. 54, 4953 (2006)]. 

PACS numbers: 61.72.Mm, 64.70. D-, 81.30.Fb, 05.40.Ca 



I. INTRODUCTION 

Grain boundaries (GBs) strongly influence the me- 
chanical behavior and other materials properties. For 
this reason, they have been widely studied both 
experimentally^ and computationally^ for decades. At 
high homologous temperature, GBs can display pro- 
nounced disorder, manifested in the most extreme case 
by the formation of nanometer-scale intergranular films 
with liquid-like properties. The formation of those films 
below the bulk melting point, typically referred to as GB 
premelting, can dramatically reduce shear resistance and 
lead to catastrophic materials failure. This phenomenon 
is of interest for predicting the formation of solidification 
defects associated with the formation of those intergran- 
ular films, which can lead to hot cracking during the late 
stages of solidificatiorPHS^ and more generally for under- 



standing the microstructure and mechanical behavior of 
structural alloys at high homologous temperature. 

GB premelting has been widely studied 
experimentally^^ as well as theoretically. Theo- 
retical approaches include discrete lattice models^ 
and molecular dynamics (MD) simulations^^, as 
well as conventional phase-field mo del j^^, wh ich 
either exploit an orientational order paramete r^ 4 * 25 ! or 
multiple phase-field d 26 * 27 ! to distinguish bet ween grains, 
and the phase-field-crystal (PFC) methocP^I], which 
resolves the crystal density field on an atomic scale and 
hence naturally models crystal defects such as isolated 
dislocations and GBs. 

Two central issues in GB premelting have been the 
determination of the premelted layer width W and the 
quantification of the fundamental forces that control this 
width. It has been generally difficult to address those 
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issues experimentally due to the challenges inherent in 
observing internal materials interfaces such as GBs. Lim- 
ited observations to date support the existence of a 
nanometer-thick premelted layer in pure materials a few 
degrees below the bulk melting point and there is more 
ample evidence for premelting in alloys. More definitive 
insights into those issues has b een provided by recent 
PFCpn anc j s i mu latior J 21 1 22 1 studies that have char- 
acterized quantitatively the structural forces underlying 
GB premelting. Quantification of those forces has been 
obtained by computing the "disjoining potential" V(W) 
defined through the excess Gibbs free-energy per unit of 
grain boundary area 

G exc (W,T) = AG(T)W + 2 lsl + V(W), (1) 

where AG = G s — Gi is the bulk Gibbs free-energy differ- 
ence between liquid (G;(T)) and solid (G S (T)) and 7,,; is 
the solid-liquid interfacial free-energy. With this defini- 
tion, y(W) represents the part of this excess due to the 
overlap of crystal density waves from the two grains on 
each side of the GB. Hence the derivative —dV(W)/dW 
measures the force between crystal-melt interfaces due 
to this overlap, which can be either repulsive or attrac- 
tive depending on whether the sign of —dV(W)/dW is 
positive or negative, respectively. It is worth noting that 
the disjoining potential can be defined in different ther- 
modynamic ensembles. While Eq. ([I]) defines it in the 
Gibbs ensemble for MD simulations, we shall also use 
the grand canonical and canonical ensembles for PFC 
and amplitude equation simulations, respectively. 

In the simplest formulation of the disjoining potential, 
W can be viewed as the width of a liquid layer sand- 
wiched between two atomically sharp solid-liquid phase 
boundaries^. Furthermore, V(W) is assumed to have a 
simple exponentially decaying form (see e.g. Ref. l30|) 

V(W) = A 7 exp(-W7<5), (2) 

where A7 = 7g b — 27^ is the difference between the GB 
energy ( 7s b) and the excess free-energy of two separated 
solid-liquid interfaces (2 7s ;). Substitution of this form 
into Eq. ([!]) and minimization of G exc (W, T) with respect 
to W predicts that, for A7 > 0, the GB is "dry" (W = 0) 
up to a bridging temperature T& < Tm beyond which 
W increases and diverges logarithmically at the melting 
temperature Tm, while for A7 < 0, the GB remains dry 
up to a superheated temperature T* > Tm', for T > T* , 
the dry GB state becomes thermodynamically unstable 

and the solid melts completely. 

Recent PFCP and MD studies! 211221 have shed light 
on several important aspects of the disjoining poten- 
tial and GB premelting that are not predicted by this 
simple model. First, while V(W) is indeed found to 
be purely repulsive for some high-energy GBs, w ith a 
form that can be approximately fitted by Eq. Q 21 29 
it is not purely attractive (i.e. dV(W)/dW > for all 
W) for lower energy GBs. For the latter GBs, V(W) 
exhibits a long-range attractive part, due to the over- 
lap of the decaying tails of crystal density waves from 



each grain, as qualitatively predicted by Eq. (J2J, but 
also a short-range repulsive part associated with the for- 
mation energy of crystal defects (e.g., dislocations for 
low-angle GBs) when W becomes comparable to the lat- 
tice spacing**!. As a result, for GBs where W does not 
diverge at Tm, V(W) generically exhibits a minimum, 
which implies that the premelted layer width remains 
finite at the melting point. Importantly, this nanometer- 
scale width represents a compromise between the long- 
range-attractive and short-range-repulsive parts of the 
disjoining potential that are both structural in nature, 
since they involve the decaying tails of the crystal density 
field and crystal-defect formation energy, respectively. 

More generally, V(W) also contains an attractive part 
due to London dispersion forces that are not accounted 
for in PFC and MD simulations, but play an important 
role in other systems such as ceramic materials^. How- 
ever, in metallic systems, dispersion forces can be es- 
timated to only contribute an attractive tail to V(W) 
whose magnitude is less than a mJ/m 2 for W on the 
nanomete r scal e. In contrast, MD computations of V(W) 
in pure N l 21 l 22 l show that V(W) has a magnitude of tens 
of mJ/m 2 for W in this same range (consistent with the 
prediction of Eq. (j2J) in a purely repulsive case). Hence, 
GB premelting appears clearly dominated by structural 
forces in metallic systems, as further supported by the 
present study in pure Fe. 

Recent PFd 2 ^ and MD^ studies have also shed light 
on the condition under which a GB will be "wet" (with 
a diverging W) or "dry" (with a finite W) at the melt- 
ing point (neglecting dispersion forces). The classic wet- 
ting condition A7 = 7g b — 27,,; > turns out to give 
a grossly inaccurate prediction when a low-temperature 
GB energy is used to compute 7^. This inaccuracy orig- 
inates from the fact that the GB energy at the melting 
point 7 3 &(Tm) is generally significantly smaller than at 
low temperature, with a large part of the decrease due to 
the elastic softening of the material (decrease of the elas- 
tic constants) at high homologous temperature 2 ^. Hence, 
a more accurate prediction of when a GB will be wet or 
dry has been obtained in PFC simulations by using the 
wetting condition A 7 > in conjunction with a value of 
lgb(TM) that takes into account this elastic softening^. 
The extension of this approach to MD simulations of Fe 
for symmetric [001] tilt GBs has yielded reasonably good 
predictions of which GBs will be fully wet or dry over a 
complete range of misorientation 2 ^. This study also asso- 
ciated departures from this prediction to the existence of 
a dislocation-pairing transition, which provides an addi- 
tional means for the GB to lower its free-energy, distinct 
from elastic softening. 

The formation of an integranular liquid-like layer 
would be expected physically to lead to dramatic de- 
crease of shear resistance of a wet GB, and hence relative 
sliding of the two grains when a shear stress is applied 
parallel to the GB plane. In contrast, a low angle dry GB 
with intact solid bridges in between isolated dislocations 
with premelted cores would be expected to support shear 
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more like a static solid. However, both theoreticaP^U 
and experimental^^! studies over the last decade have 
shown that a dry GB generically moves normal to itself 
under an applied shear stress. This motion coupled to 
shear, commonly referred to as "coupled motion" , is char- 
acterized by a relationship v\\ = j3v n between the transla- 
tion velocity of the two grains parallel to the GB plane, 
Uii, and the GB velocity normal to this plane, v n . For 
symmetric tilt GBs, analytical predictions for the cou- 
pling factor (3(0) (where 9 is the misorientation) based 
on geometrical arguments have b een va lidated by both 
MD simulations^! and experiment d 37 | 38 l 

An important question for mechanical behavior at high 
homologous temperature is when will a GB exhibit cou- 
pled or sliding motion. An MD study of symmetric tilt 
GBs in pure Np2 showed that the range of misorientation 
where sliding is observed increases close to the melting 
point, consistent with the view that the formation of an 
intergranular liquid-like film favors sliding over coupling. 
However, a recent combined PFC and MD studySS also 
showed that asymmetrical GBs can exhibit sliding unre- 
lated to premelting, so that the transition from coupling 
to sliding is generally more complex. 

In this paper, we investigate GB premelting and shear- 
ing and their relationship usingthe PFC metho d 40 " 46 !. 
We use the simplest PFC modePSELJ with the same free- 
energy function as the Swift-Hohenberg model of pat- 
tern formation 4 -^, which favors hexagonal and bec order- 
ing in two and three dimensions (2D and 3D), respec- 
tively. This model can be interpreted^ as a consider- 
ably simplified version of classical density function the- 
ory (DFT^MIH where the crystal density field is dom- 
inated by the set of primary reciprocal lattice vectors. 
With suitable choices of parameters for Fe (see Ref. HI]) , 
this model has been shown to predic t reasonable values 
of solid-liquid interfacial energ y 44 * 54 ! and GB energies^ 6 -!, 
despite the uncontrolled truncation of many other sets 
of reciprocal lattice vectors. It has also proven capable 
of reproducing a subtle dislocation-pairing GB structural 
transition at high homologous temperatures^, previously 
evidenced in a 2D PFC study of GB premelting^!, and 
also observed in MD simulations®!. However, the only 
two PFC s tudie s of GB premelting to date have remained 
qualitative^ 1321 . The study of Berry et a/.- 28 in 3D did 
not compute the disjoining potential and was not car- 
ried out for parameters of Fe. The one of Mellenthin et 
alW^ computed this potential, revealing the existence of 
a shallow minimum for dry boundaries, but was limited 
to 2D hexagonal ordering. In addition, PFC studies of 
GB shearing have only been carried out recently for 2D 
square ordering^. 

Here we compute quantitatively the disjoining poten- 
tial for [001] symmetric tilt grain boundaries over the 
complete range of misorientation < 9 < 90° and also 
study quantitatively the response to an applied shear 
stress, distinguishing between regimes of coupling and 
sliding as a function of 9 and temperature. We also 
investigate the role of thermal fluctuations on the dis- 



joining potential by carrying out PFC simulations with- 
out and with the addition of Langevin noise that is un- 
corrclatcd in space and time. Following the standard 
approach, the magnitude of this noise is fixed by the 
fluctuation-dissipation theorem and we also introduce a 
short-wavelength cut-off of this noise, which is shown to 
be necessary to avoid unphysical divergences. To bench- 
mark our results, we compare the disjoining potentials 
computed by PFC and MD simulations. The MD simu- 
lations are carried out using the same EAM potential 
that was previ ously us ed to calibrate PFC model pa- 
rameters for F (j 44 l 45 l 54 J With appropriate rescaling of 
length and energy, the PFC model can be expressed in 
a form that involves a single small dimensionless pa- 
rameter e. In the context where PFC is derived from 
DFT, e is u niquely determined by liquid-structure fac- 
tor propertie d 44 * 45 * 54 !. More generally, decreasing this pa- 
rameter makes the freezing transition more weakly first 
order and increases the width of the solid-liquid inter- 
face in units of lattice spacing. To gain additional in- 
sights into the role of this parameter, we also compute 
the disjoining potential using amplitude equationa^HSSl. 
These equations can be formally derived from the PFC 
model in the small e limit using similar multiscale expan- 
sions introduced previo usly to analyze continuum models 
of pattern formatiorP^HHl While amplitude equations 
have recently been shown to break down for high an- 
gle GBs because of issues related to frame invariance^!, 
they are asymptotically exact for small e and low angle 
GBs. Hence, in the present context, they allow us to 
derive scaling laws for the range of misorientation over 
which GBs exhibit a diverging liquid-layer width. Here 
we report amplitude-equation results that pertain to the 
numerical computation of disjoining potentials and the 
derivation of this scaling law. In a companion paper to 
this one, we present the results of a more in depth an- 
alytical and numerical amplitude-equation study of dis- 
joining potentials. In this paper, we also treat the case 
(analogous to the computation of 7 surfaces) where two 
crystals of the same orientation are translated from each 
other in the plane of the interface^. 

The rest of this paper is organized as follows. In sec- 
tion [Hj we write down the PFC model and motivate the 
necessity of introducing a short- wavelength cut-off in the 
Langevin noise with a short calculation of the renormal- 
ization of the free-energy of the liquid state. In section 
III we write down the amplitude equations and relations 
that connect them to the PFC model. In section [rv] we 
then discuss the methodology to compute the disjoining 
potential with the PFC model without and with fluctu- 
ations and in MD simulations. The same methodology 
developed in Refs. [21] and 122] which relates the disjoin- 
ing potential V(W) to the probability density P{W, T) of 
width fluctuations at different temperatures close to the 
melting point is used for both MD and PFC simulations 
with noise. PFC simulations without noise follow the 
same methodology developed in Ref. 29 using the grand 
canonical ensemble appropriate for PFC. In section[V] we 
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then present the numerical results for disjoining poten- 
tials computed with the PFC approach without and with 
fluctuations, MD simulations, and amplitude equations, 
and we derive a scaling law for the critical angle below 



(or above) which a GB is dry or wet. Next, in section VI 



we present the PFC results for GB shearing. Our results 
are then summarized in section IVHl 



II. PFC MODEL 



A. Basic equations 



We use the simplest version of the PFC modePHHwith 
the same free-energy functional as the Swift-Hohenberg 
model of pattern formation^, which is defined by 



I = / df 



2\2i 



: [a + X(q 2 +V 2 ) 



(3) 



where <f> is a dimensionless measure of the crystal density 
field. This free-energy favors bcc ordering in 3D and has 
been shown to predict reasonably well the solid-liquid 
interfacial free-energy and its anisotropy as compared to 
MD simulations for parameters of Fe^H 

In order to minimize the number of relevant parame- 
ters, we introduce the scalings 



(4) 

(5) 
(6) 

(7) 



if, 



F 




leading to the dimensionless free-energy functional, 



F = df 



(8) 



where we have dropped the prime onr', so that r denotes 
the dimensionless position vector hereafter. 
We use a non-conserved dynamics of the form 



dtp 



sn 



(9) 



where Q = F — fi J drip is the grand potential, fi is the 
constant chemical potential, and r\ is a Langevin noise 
that is uncorrelated in space and time with zero mean 



(10) 



and variance 

Mf.tW.O) =rS(r-r<)S(t-ti), (11) 

where 

2k B Tg 



r = 



(12) 



measures the noise strength in our dimensionless PFC 
units. In simulations without noise, 77 is set to 0. 

While the present study could be equivalently carried 
out using t he st andard form of the PFC dynamics that 
conserves we have found it useful to work in the 

grand canonical ensemble where [i plays an analogous 
role to temperature in the Gibbs ensemble in which the 
disjoining potential is usually defined via Eq. ([I]). We 
can then use thermodynamic relations to relate results 
in those two ensembles. 

The higher order spatial derivatives make the PFC 
equations stiff in real space. We therefore conduct our 
simulations in spectral space, as originally presented irP^ 
with the details for the addition of noise presented in 
Appendix |A 1| 

To simulate GB shearing, we stay in the grand canon- 
ical ensemble, but use the modified PFC model with in- 
ertia (i.e. with an additional d 2 ip/dt 2 term)P^ 



o 



<9V 

at 



5F 

'Jtp 



(13) 



which has the advantage to quickly relax the elastic field 
via propagative phonon-like modes, as opposed to purely 
diffusive modes. This allows us to shear the crystal at 
velocities much faster than those we could use with the 



standard diffusive dynamics. Eq. (13) is simulated with 
a modified spectral scheme presented in Appendix |A 2| 



B. Short-wavelength noise cutoff 

In this section, we examine how to choose the short 
wavelength (high k) cutoff of the Langevin noise in the 
PFC model. Previous studies of Langevin noise in PFC 
have used such a cutoff to carry out simulations^^. 
The issue of whether introducing noise may lead to a 
double counting of fluctuations has often been discussed. 
On the one hand, the free-energy functional of classi- 
cal DFT, or its simplified version in PFC, already rep- 
resents a mean-field average of atomic scale fluctuations. 
Thus adding back fluctuations would appear to be dou- 
ble counting. However, without the addition of noise, 
dynamical DFT or PFC simulations do not reproduce 
long-wavelength hydrodynamic fluctuations, e.g. capil- 
lary fluctuations of the solid-liquid interface^!, which are 
present in a real system. Hence, in situations where those 
fluctuations, which are not double counted, are of inter- 
est, it makes sense to add noise with a short- wavelength 
cutoff to avoid double counting short-wavelength fluc- 
tuations already incorporated in the DFT or PFC free- 
energy functional. 

In general, the magnitude of the noise is set by the 
fluctuation-dissipation theorem. However, the choice of 
the cutoff between not-double-counted long-wavelength 
and double-counted short-wavelength fluctuations con- 
tains some degree of arbitrariness. A common choice is 
to cut off the noise on the lattice scale, which is equiva- 
lent to choosing the maximum wavector of the noise k max 
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of order unity in our dimensionless units where length is 
measured in unit of q^ 1 . To make this choice somewhat 
more precise, we compute here explicitly how the noise 
renormalizes the grand potential (Landau free-energy) of 
the liquid. This renormalized noise-induced excess po- 
tential is found to diverge ~ k^^. We then compute the 
value of k max at which this excess potential is equal to the 
barrier height of the double well potential between solid 
and liquid of the noiseless PFC. This value sets an upper 
bound on k max , beyond which the noise is too strong and 
essentially destroys solid-liquid coexistence. 

To compute the renormalized grand potential, we sub- 
stitute rp — $i+Stp into Eq. Q and expand to quadratic 
order the excess An = n^^ — nfyi), which is the differ- 
ence between the grand potential of a liquid with density 
fluctuations driven by the Langevin noise and a liquid of 
constant density tpi. Since \i J dr5%p = gives no contri- 
bution, this excess is given by 



Aft = df 



(-e + 34>f + (V 2 + l) 2 ) 5^{r) 



dk 



2(2tt) : 



7 j(k)Sip(k)Sip(-k), 



where we have defined 7 (A;) 
the Fourier transforms 



(14) 

+ 3$ + (1 - fc 2 ) 2 and 



6${k) = 
Stjj(r) = 



drStp(r)e- lk - F , 
1 



(2tt) e 



ik-r 



dk8ij){ky K - 



(15) 
(16) 



Next, the "renormalized" excess grand potential is the 
equilibrium statistical average of An given by 

/dk ~ - - 

_3 7 (fc)(^(A;)^(-*))e 9 . (17) 

The correlation function lM>{k)5^i{—k)) eq is readily ob- 
tained by linearizing Eq. (|9| around the liquid state with 
the substitution ip(r, t) — ipi + 5ip(f, t) and Fourier trans- 
forming the resulting equation, which yields 



dtS^(k) = -j(k)8ip(k) + fj(k, t). 



Eq. ( 18 ) has the solution 
rt 



Si>(k,t)= dre~ 7 W(*- T )?7(fc,T). 



(18) 



(19) 



Substituting this form into Eq. ( 17 ) and using the Fourier 
transform of Eq. ( 11 ) 



(rj(k,t)f)'(k',t')) =r(2Tr) 3 S(k + k')5(t-t'), (20) 
we obtain 

(27r) 3 r 



(Sip(k,t)6ij}(k' ,t)) 



2 7 (fc) 



(l-e- 2 ^) S(k + k'), 



(21) 



and thus in the equilibrium long time limit (t — > 00) 



(22) 



Substituting the above result into Eq. (17 1, we obtain 



the final form of the excess liquid grand potential due to 
noise-induced density fluctuations 



(An) t 



4nTV 
4(2tt) 3 



k 2 dk 



4irTV 3 
12(2tt) 3 max ' 



(23) 



Eq. ( 23 1 shows that the renormalized grand potential 



of the liquid has an ultraviolet divergence that requires 
keeping k max finite. Next, to estimate what this cutoff 
should be, we compare (AQ) eq to the barrier height of 
the double-well potential between solid and liquid. To 
compute the latter, we start from the expression 



AQ = Q, 



fii 



where 



W=(-e + l)f + f 



(24) 



(25) 



and n s /V is obtained by substituting into Eq. ([8| the 
form for a bec crystal density field 4 ^ 



ip(r) -tp = e 1/2 ^2 Ae lki ' F = 4e 1/2 A(cos qx cos 



<l 1 J 



■ cos qx cos qz + cos qy cos qz) 



(26) 



where fej are the set of 12 [110] bec principal reciprocal 
lattice vectors (cf, Eq. ( 33 1 below) and q = 1/ y/2. This 
yields the expression 

n s /v = (-e + i) y -f + ^--^ s 

-6e 2 A 2 + 18eip 2 A 2 + 48e 3/2 ?M 3 + 135e 2 A 4 . (27) 

Expanding tp s ^ in powers of e 1 / 2 as ip s ^) = ip s (i)oe^ 2 + 
ip s (i)i e + • • • and using the results (see Section III A of 

Ref.SH) that $8(1)0 = V>c = -^45/103 and ^ s (/)i = 0, we 
obtain the expression for the grand potential difference 
between solid and liquid per unit volume 



An = e 2 Vf dw (A)+0(e 3 ), 



(28) 



where 



fdw(A) = ~6A 2 + 18iJj 2 A 2 + 48ip c A 3 + 135 A 4 (29) 

has the shape of a double well potential with minima at 
Ai = and A s = 8/ (3V515) corresponding to the liquid 
and solid states, respectively, and with a maximum at 
Ab = 4/ (3v515) corresponding to the free-energy barrier 
between those two states. Next, equating the renormal- 
ized excess grand potential of the liquid and the height 
of the double-well potential, (An) eq = e 2 Vfd w (Ab), and 



() 



using Eqs. (12) and (23) together with the expressions 



for g, A, and e given in Appendix C, we obtain 



,3 _ 3(2^) 2 iV-ff( go )C"( go )3^n 
where we have defined the dimensionless constant 



B 



(31) 



Finally, evaluating fdw(Ab) using Eq. (29) and using 



the parameters given in Table [XTJ we find that k max m 
0.918, which implies that the shortest wavelength of the 
noise converted back to dimensional units is \ m in = 
2tt / ((Journal) = 0.77a where a is the unit cell size. Cut- 
ting off the noise at this lengthscale ensures that the noise 
does not wash out the barrier between the solid and liq- 
uid states. To be safe, we restrict ourselves to cutoffs 
on the order of one unit cell or larger; details are pre- 
sented in Appendix |A 1| This choice clearly affects the 
relative impact of the noise, and we present results for 
two different choices of X m in below. 



III. AMPLITUDE EQUATIONS 

As a complementary approach we use here an ampli- 
tude equations (AE) model which is derived via a mul- 
tiscale expansion from the three-dimensional phasc-ficld- 
crystal model, Eq. ([8]), and use e as a small parameter 
in the regime of the phase diagram which describes the 
coexistence between the bcc and the homogeneous (melt) 
phase. 

As in classical density-functional theory (DFT), the 
spatial variations of the density field, 5ip(r), are expanded 
as a sum of density waves 



(j) ik ij >-f 



(32) 



where ki represent different reciprocal lattice vectors 
(RLVs) and v/^ are their associated amplitudes. In the 
liquid phase, where the time averaged density is spa- 
tially constant, the amplitudes vanish, and in an undis- 
torted solid phase they all attain the same constant value, 
u 1 -^ = u s . 

Unlike DFT, where a large number of modes is required 
to obtain sharp peaks around atomic positions, the sim- 
pler free-energy allows for a truncation of this sum to a 
small set of reciprocal lattice vectors. Various methods 
have been developed to change the kernel of the free- 
energy in order to stabilize a variety of two and three di- 
mensional periodic and crystal structures. Here we focus 
on bcc structures, therefore we restrict the summation to 
the N — 12 principal reciprocal vectors 



[110], [101], [Oil], [110], [101], [Oil] 
[110], [101], [011], [110], [101], [011]. 



(33) 



Notice that by the condition of having a real density field 
ij) the N complex amplitudes are not independent but 
are complex conjugate (denoted by a star) for antipar- 
allel RLVs. Therefore, we restrict the description to the 
first row of these RLVs and use only N/2 independent 
complex fields vP\ 

A detailed derivation of the amplitude equations, 
which describe the evolution of the fields has been 
given in Refs. EH and 123 and therefore we only give the 
resulting expressions here. The evolution equations can 
be derived from a free-energy functional, which reads 



l-JV/2 N/2 

Fae =F dR J2 Pi^ W | 2 + To E A{j)Aii> 
■I -' - J= i 




\ 2 N/2 



'N/2 

+2Al w A*i AiQiA w i + 2Ano^4iio^4*oi^ioi 
+2^110^011^011^110 + 2 ^iio^oii^oii^iio 



-2A 011 A^ oI ^4ioi^4oii + 2A* 011 A W1 A* 101 A 011 j 

■g (^011^101^*10 + Aui^-ioi^ilo + -^011^110^*01 

-A 011 Al w A w i + A* 01l Ai W A* 101 + A 01 iA* 110 A 101 



"^011^101^110 + -^011-^101-^11 



Ft . 



(34) 



Here, we have introduced rescaled amplitudes 

A (j) = n (3)/ Us (35) 

and a dimensionless "slow" scale R which is introduced 
below. 

In a more general context the above expression can be 
derived from classical DFT and yields 



Fae 



n Q k B T 



dr 



(X 



t U) u U> 
S(q ) 



C"(q ) 



(36) 



where the function f({u^}, T) contains the higher order 
nonlinear terms in the amplitudes u^' and an explicit 
dependence on the temperature T . The form of these 
nonlinearities depends slightly on the underlying model: 
Above it is given for a PFC model, and there are some 
differences if these terms are derived from DFT using an 
equal weight ansatz. However, the differences are small 
and lead e.g. only to tiny changes of the anisotropy of 
the surface energy, as had been investigated irP^. Fur- 
thermore, no is the density in the liquid state, C(r) is 
the direct correlation function with Fourier transform 



(37) 



C{q) = no drC{r) exp(— iK ■ r) 
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S(q) = 



(38) 



with r = \r\ and q = \K\; it is related to the liquid 
structure factor by 

1 

Here, all quantities are evaluated at the (first) peak of 
the structure factor qq. We introduce a (dimensionless) 
length scale 

R = e 1/2 q f 

with 

24 96 



S(q )C"(q )q 2 103 
This allows the identity 



(39) 
(40) 



F 



noksT „ _ x 2 — 1/2 
C (q )q a u s e ' . 



The differential operator \3j is given by 



-v 2 , 



(41) 



(42) 



where the nabla operator acts on the slow scale R, and 
the vectors few) are the normalized principal RLVs. The 
second term in the operator preserves the rotational in- 
variance of the equations and is related to the use of the 
nonlinear strain tensor in elasticity. 
The thermal tilt 



Fn 



T-T, 



^ £-3/2^-3 



j d Rj: 



(43) 



is added phenomenologically to the model. Here, Tm is 
the melting temperature, L the latent heat and <p is an 
"order parameter" which discriminates between solid and 
liquid, with 



0({A»}) = |^v / AWA( 



N/2 



j)* 



(44) 



We note that this expression is invariant under elastic 
deformations and lattice rotations, which affect the com- 
plex phases of the amplitudes. The choice of the coupling 
function and its implications are discussed in more detail 
in Ref . [7H 

Thermodynamic equilibrium corresponds to a station- 
ary state of the free-energy functional, and we use relax- 
ation dynamics according to 



(45) 



dA® _ 6F 
dt ~ " j 6lU>' 
Since we focus here exclusively on static situations, the 
choice of the kinetic coefficients Kj is arbitrary. 

This description predicts the correct anisotropy of 
surfaces energies^ and elastic properties and contains 
naturally the linear theory of elasticity. The above 
Ginzburg-Landau model for crystals is conceptually close 
to theories of superconductivity and pattern formation in 
hydrodynamics^. 



IV. GRAIN BOUNDARY PREMELTING: 
METHODOLOGY 

A. Disjoining Potentials without noise 

While the definition of the disjoining potential in the 
introduction of this paper is framed in terms of a Gibbs 
ensemble (constant T, p, and N), it is easier in the PFC 
model to conduct simulations in the grand canonical en- 
semble (constant T, /i, and V). As is discussed in Ref. [2^1 
studying grain boundary behavior as a function of fi—fieq 
is equivalent to using T — Tjj. In the spirit of Eq. 
we can define the disjoining potential as the excess grand 
potential energy in the system 

= L y L z [(L x ~W^))uj s ^)+W{^ l (p)+2 lsl +V(W)] 

(46) 

where 0J s (i) (/u) are the grand potential densities of of the 
solid (liquid) respectively, W is the width of the liquid 
layer, L V L Z is the area of the grain boundary and L x is 
the length of the system perpendicular to the boundary. 
We can then define the grain boundary energy 



fi(/i) - L x L y L z u s (n) 



LyL z 



-2 lsl + V[W(»)}, 



(^(m)-^(m))^(m) 

(47) 



In order to calculate the liquid layer width W we relate 
it to the excess mass in the system due to the presence 
of a grain boundary. We define an excess mass per unit 
area of the grain boundary by comparing the mass of two 
systems at the same chemical potential, one a perfect 
solid and one containing a grain boundary, 



(48) 



By also obtaining the liquid density at the same chemical 
potential we can say that the excess mass must be equal 
to the liquid thickness times the difference in densities of 
the liquid and solid phases 



W[i)i((jL) - ip a (p)] = Vw(m)- 



(49) 



Combining the two equations gives a definition for the 
liquid layer W(p), 



W(jj) = L, 



(50) 



In order to calculate the disjoining potential from 
Eq. (47), we need a way to calculate the grain bound- 
ary energy 7o b . By dividing Eq. (j46j> by L x L y L z and 



using Eq. (47) we see that 



u = oj s (fj,) +7 gb (/i)/L 2 



(51) 



As seen in Fig. [T]we can then obtain the grain bound- 
ary energy by running a series of simulations at fixed 
chemical potential and varying system size. Plotting to as 
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FIG. 1. Grand potential energy density ui at a constant chem- 
ical potential p = —.19609 for three different inverse system 
sizes 1/L X perpendicular to the grain boundary for a mis- 
orientation of 9 — 18.9°. The intercept is twice the grain 
boundary energy due to periodicity. 



a function of the inverse system size, we see a clear linear 
relationship, and are able to extract 7„5. We should also 
note that as the system is periodic we have two identical 
grain boundaries and the resulting energy should be di- 
vided by two. Theoretically, all that is needed to compute 
the disjoining potential is to calculate the grain bound- 
ary energy for a whole range of chemical potentials, and 
then use Eq. (47) to extract the disjoining potential for 
each fi. However, in practice this is computationally un- 
feasible. Instead we choose an alternate route and define 



the disjoining pressure II 

i on 



L y L z dW 



V'(W). 



(52) 



We also recognize that the grand potential densities ui 
can be expanded as a function of [i, 



uj s - u)i « -(ips 9 - ipi q ){n ~ Me?)- 



(53) 



When the system is at equilibrium II is equal to zero and 



(54) 



As noted in Ref. 29, V'(ji) is a known function of /i which 
depends only on bulk thermodynamics. By conducting 
simulations for a range of fixed chemical potentials we 
can obtain V'(W) which is then integrated to get V(W) 
as long as we know 7 g (, for one particular fi (calculated 
from Eq. (51 1). This method is far more computationally 
tractable than calculating V(W) directly from Eq. p7| . 

While this discussion as focused on undercooling as a 
function of chemical potential (i, it is possible to relate 
this quantity to undercooling as a function of tempera- 
ture. As derived in the supplementary material of Ref. l23l 
we relate the energy difference between solid and liquid 
phases to the difference in grand potential energy densi- 
ties and expand oj in fx 



L- 



T -T 



M 



T, 



UJ S 



-OA 



e? 



(55) 



where the latent heat of fusion for the MD Fe potential is 
L = 1.968 • KT 24 kJ/A 3 . This allows us to present PFC 
results as a function of temperature rather than chemical 
potential. 

For the amplitude equations we use a similar approach 
to define the melt layer thickness and the disjoining po- 
tential. Here we focus on the case that the thermal tilt in 
the free energy is a linear function in amplitudes, which 
has the consequence that the bulk amplitudes in solid 
and liquid deviate from and 1, respectively, which are 
the potential minima for T = Tm- In the bulk all am- 
plitudes have the same magnitude A (they differ only in 
phase due to the grain rotation), and therefore the free 
energy density becomes 



f AE (A,T) = f (A) + f T 

= Fo? /2 q 3 o ( U 2 - A 3 + l 



A 4 



AL 



T-T M 
T, 



Notice that in the bulk the gradient square term does 
not appear despite the spatial variations due of the com- 
plex phases; this is a property of the box operator, as 
discussed in detail in Ref. 63] The minimum positions of 
the bulk free energy density can be found as the roots 
of the cubic equation f' {A) + f T {A) = 0, and we denote 
them by Ai(T) and A S (T) (the third root is the maxi- 
mum in between these two minima). We define the melt 
layer thickness for a two-dimensional system as 



NL„L, 



3 = 1 



where the interface area L y L z is measured on the slow di- 
mensionless scale. The prefactor g^e -1 / 2 stems from the 
conversion of W to the dimensional "short" scale. We let 
the simulations relax towards equilibrium with the bias 
potential Ft, which does not enter into the expression of 
the disjoining potential y(W). The latter is calculated 



V(W) = 



LyL z 



dRf - LyLzWfiAt 



L y L z (L x - W)f{A s 



(57) 



in analogy to Eq. ( 46 ) above. 



By tuning the bias potential Ft we can therefore ex- 
tract the convex parts of the disjoining potential, as be- 
fore for the PFC simulations; this approach is only possi- 
ble in the convex regions, V"(W) > 0, otherwise the equi- 
librium state is unstable. To also get estimates for the 
concave parts for the attractive potentials we instead do 
dynamical runs at T = Tm, starting with a rather wide 
liquid layer. Due to the attractive interaction the melt 
starts to solidify, and we measure W and F as function 
of time. Although the system is not in full equilibrium 
during this evolution, the extracted disjoining potential 
V(W) matches well to the static results which were ob- 
tained by the above method. 
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B. Thermal fluctuations 

Calculating the disjoining potential from a fluctuat- 
ing interface requires a different technique then presented 
above for PFC with out n oise. Instead we follow the ap- 
proach of Hoyt et al! 21 * 22 *, which relates the excess inter- 
facial free-energy to fluctuations of the liquid layer width 
as a function of undercooling. The grain boundary free- 
energy 7 g f,(/i) is defined as above in Eq. (17), allowing us 
to easily write the probability of observing a liquid layer 
width W as 



P(W,T) = Cex.p[-A lgb (T)/k B T], 



(58) 



where A is the interface area, C is a temperature depen- 
dent normalization constant, and we have used Eq. (55) 



to convert from fi to T. We obtain these probability dis- 
tributions by conducting simulations with noise for vari- 
ous undercoolings and plotting the resulting histograms, 
as seen in Fig. |16[ The inclusion of noise prevents us 
from using the method presented above to calculate the 
liquid layer width, due to the large fluctuations in aver- 
age density tj). Therefore we use the procedure presented 
below in Appendix [B] to define the width W . Once we 
have the width histograms we can extract the disjoining 
potential by writing 

V(W) = -(k B T/A) \nP(W,Ti) - Au> f W - 2 lsl + a t , 

(59) 

where i enumerates simulations run at different temper- 
atures Ti, Aujf is the bulk grand potential energy den- 
sity difference between solid and liquid for a particular 
undercooling, and the ai are free constants. The dj are 
then used to match the various disjoining potential curves 
(from different undercoolings) into one continuous dis- 
joining potential as seen in Fig. [17] 



C. Molecular Dynamics simulations 

For the purpose of comparing with PFC and AE pre- 
dictions, MD calculations of disjoining potentials have 
been undertaken for bec Fe using the method described 
in Refs. [2"T1 and [2"2"1 This approach for computing 
by MD involves the use of equilibrium MD simulations 
to compute the grain-boundary width histograms, from 
which the disjoining potential can be extracted using 



equations analogous to Eqs. (58) and (59) (see Ref. 1211 for 
details). In the current application of this approach we 
performed MD simulations using the LAMMPS codtP. 
The simulations were performed in the NP Z AT ensem- 
ble, corresponding to fixed particle number N and cross- 
sectional area A, with the temperature T and pressure 
normal to the grain boundary P z — controlled by a 
Nose-Hoover thermostat and barostatP^. Use was made 
of periodic boundary conditions, such that the simulation 
system contained two identical symmetry tilt boundaries. 
Interatomic interactions were modeled by the embedded- 
atom-method potential given in Ref. [5H1 Width his- 



tograms were computed from snapshots of equilibrium 
MD simulations lasting at least 80 ns, at two or more 
temperatures below, but within 10 K, of the melting 
point. For each snapshot, the grain-boundary width was 
determined through the use of a structural order param- 
eter based on a formulation due to Morris^. Further 
details can be found in the supplementary information 
to Ref. H3 



V. GRAIN BOUNDARY PREMELTING: 
RESULTS 

In this section we discuss the results for grain boundary 
premelting in the context of phase field crystal (PFC) 
and amplitude equations (AE) modeling. In particular, 
we inspect symmetric tilt grain boundaries in bec S iron 
for misoriented [100] surfaces. The results are compared 
to MD simulations. 

There are two relevant sets of information that are ex- 
tracted here. 

First, the width of the melt layer W forming at the 
grain boundary as function of temperature is discussed. 
In thermodynamic equilibrium, an infinitely wide liquid 
phase can coexist with the two grains, and in this limit in- 
teractions between the grains with different orientation 
have decayed completely. Therefore, the curves W(T) 
must diverge at T = Tm- For a repulsive grain-grain 
interaction this leads to monotonously increasing and di- 
verging functions W(T) from below the melting point. 
For attractive grain boundaries, the dependence is not 
unique, since the diverging branch is located above the 
melting point and unstable. Additionally, a stable branch 
with finite melt layer thickness at T = Tm exists with 
lower W, and it merges with the unstable branch at a 
finite overheating temperature T D > Tm- Above the 
melting point this branch is only metastable, as from 
bulk thermodynamics the system would reduce its total 
free energy by complete melting. Here, we only track 
the (meta-)stable branches of the solution, as they can 
be obtained by simple gradient dynamics. For a more 
extended discussion and the consideration of unstable 
branches in the framework of phase field modeling we 
refer to Ref. [UJ In contrast to this kind of modeling, 
which is characterized by a constant order parameter in 
each bulk phase, we preserve in the PFC and AE simula- 
tions the atomic structure, and therefore the interaction 
can be studied as function of the grain misorientation, 
and one finds a transition from attractive situations for 
low angle grain boundaries to repulsive interactions for 
high angle boundaries. 

Second, related to the W(T) curves is the disjoining 
potential V(W), which describes the interaction between 
the grains which are separated by the melt layer of thick- 
ness W. The disjoining potential is normalized such that 
it decays to zero for W — > oo, when the excess energy at 
the melting point reduces to twice the solid-liquid inter- 
facial energy per unit area, 7^. A repulsive interaction 
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corresponds to V'(W) < 0, and an attractive case to 
V'(W) > 0. Typically, an attractive situation at large 
distances still comes with a repulsive interaction if the 
interfaces arc very close to each other. 



In subsection V A we present the deterministic results 
of the continuum simulations, i.e. without thermal fluctu- 
ations. The deviations from the MD results, which show 
a transition from attraction to repulsion at lower angles, 
is interpreted as a different ratio of the solid-liquid in- 
terfacial energy j s i and the grain boundary energy 7 g ^. 
Therefore, in subsection |V B| this ratio is set to a similar 
value in the continuum theories, leading to a shift of the 
curves in the right direction. 

Finally, in subsection |V C| the influence of thermal 
noise is studied in the PFC model, which also leads to 
a shift of the disjoining potential closer to the atomistic 
results. 



A. AE, PFC without noise and MD results 

The PFC results for liquid layer width as a function of 
undercooling are presented in Fig [5] As detailed above, 
for each misorientation the system is relaxed at constant 
chemical potential [i and the liquid layer width is mea- 
sured. Then fi is increased and the system is relaxed 
again, until the amount of liquid diverges. PFC simula- 
tions are conducted with At = 0.5, and Ax w a/8 where 
a is one lattice spacing. In order to ensure periodic simu- 
lation geometry for rotated bi-crystals, we follow the pro- 
cedure set forth in Appendix B of Ref. [29] The system 
size in the direction normal to the boundary (L x ) is kept 
as close to L x « 100a as possible given the constraints 
of periodicity, while the system size along the boundary 
is set to the minimum length allowed under periodicity 
and L z = 2a. 

There are qualitative similarities and differences to 
sharp interface theory. Here we find that, for the repul- 
sive boundaries, there is indeed a logarithmic divergence 
of the melt layer thickness if the temperature reaches 
Tm from below, which indicates short range interactions. 
For low temperatures, there is no bridging temperature. 
Instead there is always a small liquid layer in existence. 

The disjoining potentials as integrated from Eq. (54) 
for e = 0.0923 are plotted in Fig[3j Here we can easily see 



the transition between repulsive and attractive disjoining 
potentials at 6 rj 37°. 

The corresponding results for AE are shown in Figs. [4] 
and[5j Here we typically use a lattice spacing of Ax ss 0.5 
(measured on the dimensionless slow scale) in a spectral 
code. They show a qualitatively similar behavior, and 
in particular the results of this coarse-grained theory de- 
liver results which are on the same order of magnitude. 
The main difference is that the behavior is already sig- 
nificantly more repulsive, and the transition between at- 
traction and repulsion occurs around 9 = 25°. We also 
show here results of the concave parts of the disjoining 
potential curves, V"(W) < 0, which were obtained by 
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FIG. 2. (Color online) Grain boundary liquid layer width as 
a function of undercooling for various misorientations at a 
value of e = 0.0923. There is a transition from attractive to 
repulsive grain boundaries at approximately 9 — 37° . Data is 
from PFC simulations. 
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FIG. 3. (Color online) Disjoining potential as a function of 
grain boundary width for various misorientations and e = 
0.0923 from PFC simulations. There is a transition from at- 
tractive to repulsive disjoining potentials at approximately 
6» = 37°. 
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FIG. 4. (Color online) Width of the melt layer as function of 
the deviation from the melting temperature for e = 0.0923, 
as obtained from AE. 
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FIG. 5. (Color online) Disjoining potential for inclined (100) 
interfaces for e = 0.0923, obtained from AE. Below a criti- 
cal misorientation 8 C ~ 25° the interfaces attract each other, 
above this value they are repulsive. The thin (concave) parts 
of the attractive potentials are obtained by dynamical simu- 
lations. 
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FIG. 7. (Color online) A comparison of disjoining potentials 
between PFC, AE and MD for a misorientation of 9 « 44°, 
with e = 0.0923 and e = 0.0329. The MD results represent 
statistics generated from three different undercoolings. 

as AE is derived from the PFC model by a multiscale ex- 
pansion, which assumes that the length scale over which 
the amplitudes vary spatially is large compared to the 
lattice unit. This assumption breaks down for a high 
angle grain boundary, which introduces a short-scale os- 
cillation of the amplitudes, and therefore discrepancies 
between the models have to be expected, as discussed in 
detail in Ref. 

A comparison with the MD data in Fig. [7] shows that 
the two continuum theory curves are very close together 
for a high angle misorientation of about 45°, but still 
quantitatively different from the atomistic results. The 
reasons for this disagreement will be elucidated in more 
detail in the next section. 



FIG. 6. (Color online) Disjoining potential comparison for 
e = 0.0923 for selected misorientations for PFC and AE data. 



using dynamical simulations. As mentioned before, dy- 
namic and static simulations deliver comparable results 
in the range where both of them apply; the small dis- 
crepancy can be seen best for the 11.4° misorientation, 
where only a small mismatch of the order of 1 mJ/m 2 ex- 
ists. Similarly to the stronger repulsion, the overheating 
range, i.e. the temperatures up to which a bicrystal is 
still thcrmodynamically metastable, is lower accordingly. 
A detailed comparison of the disjoining potentials, which 
were computed by these two continuum models is shown 
in Fig. [6] It turns out that the positions of the energetic 
minima are very similar, although the definition of W is 
not the same in the two models, and therefore horizontal 
shifts of the curves would be conceivable. Nevertheless, 
the energy values, in particular in the minima, are still 
quite different for low angle misorientations. 

We point out that this disagreement is not surprising, 



B. Rescaling of e 

As seen in Fig. [7J although qualitatively similar, there 
is a large difference in the range of interaction between 
MD and PFC. 

Here we aim at a better matching of the continuum 
and the atomistic simulations. To this end, we first in- 
vestigate the critical angle for the transition between at- 
traction and repulsion as function of the only free pa- 
rameter, which appears in the (rescaled) PFC and AE 
models, namely e. This dependence is shown in Fig. [HJ 

The observation is that for low values of e the criti- 
cal angle scales as 9 C ~ e 1 / 2 both for PFC and AE, and 
only for higher angles the results differ from each other. 
In particular, there is no premelting transition in PFC 
for e 1 / 2 > 0.3, in agreement with findings in Ref. H251 On 
the other hand, the critical misorientation curves are two 
separate straight lines, coming from (100) and (110) in- 
terfaces; this is a consequence of the incorrect treatment 
of the (discrete) rotational invariance of the model, as 
discussed in detail in Ref. IT21 

The basic physical argument for the premelting transi- 
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FIG. 8. Plot of critical angle 8 C as a function of e 1 ^ 2 for both 
PFC and AE. Disjoining potentials are repulsive (attractive) 
below (above) those lines. The range of misorientation be- 
tween triangles where GBs are observed to form a liquid-like 
layer in MD simulations^ is in good quantitative agreement 
with the range predicted by PFC and AE for the value of e 
that matches the ratio of bulk modulus and solid-liquid inter- 
facial free-energy in MD simulations. The AE theory predicts 
the asymptotic form 6 C ~ e 1//2 (also -n — 9 C ~ e 1//2 for misori- 
entations close to 90°) for e<l. 



comparison. For low angle grain boundaries, the spac- 
ing between neighboring dislocations scales as 1/qod. On 
the other hand, a dislocation is a defect, which corre- 
sponds to a singularity in the continuum theories. For 
a Burger's circuit around the dislocation line the phase 
of the complex amplitudes increases by a multiple of 27r, 
and at the same time density and the magnitude of the 
amplitudes tends to zero inside the dislocation (or at least 
they become small). Therefore, a dislocation core is sim- 
ilar to a solid-liquid interface, and its diameter is there- 
fore proportional to the extent of a solid-liquid interface, 
i.e. qQ 1 e~ 1 / 2 . The transition between attraction and re- 
pulsion happens when the dislocation spacing becomes 
comparable to the dislocation core size, thus 9 C ~ e 1 / 2 . 

In AE, which is derived from the PFC model in the 
limit of small values of e, the scaling 9 C ~ e 1 / 2 can also 
be inspected from another perspective. For small e the 
correction term in the box operator becomes negligible, 
and therefore e drops out of the dimensionless bulk equa- 
tions. It only enters into the problem by the geometrical 
setup of the grain boundary, which requires the rotation 
of the crystals. This leads to amplitudes 



(65) 



tion is that 2^/ s i — "f g b. The scaling of the surface energy 
is 



7s; ~ n k B T 



-C"(qo) 
%o) 



1/2 



(60) 



The grain boundary energy, on the other hand, stems 
from the elastic energy of the geometrically necessary 
dislocations at the grain boundary. For fixed misorienta- 
tion the strain is therefore fixed, and the elastic energy 
scales with the elastic constants. Thus we obtain 



B~n k B T(-C"(q ))qlu 2 s , 



(61) 



with B being the bulk modulus. The ratio of these two 
parameters therefore scales as 



B 

Isl 



9oe 



-1/2 



(62) 



For low misorientations the grain boundary energy de- 
pends on as given by the Read-Shockley law, 



lab = 1q0{A -In 6) ~ 7o 6 



(63) 



where 70 and A depend on the elastic constants, the 
Burgers vector and the dislocation core radius; the second 
relation is the asymptotic behavior for small misorienta- 
tions. Comparing this to Eq. (62) therefore implies for 
the critical misorientation with 7 s j ~ 2j s i 



e 1/2 for e < 1. 



(64) 



A more geometrical interpretation of the transition be- 
tween attraction and repulsion stems from a length scale 



on the slow scale R, with the unity matrix 1 and the 
usual rotation matrix 



R 



cos <p 
=F sin i 



± sin< 
costb 



(66) 



with (j> = 9/2. For low angles the argument of the expo- 
nential becomes 







±0 



,1/2 



,1/2 







R. 



(67) 



Hence the misorientation and e enter in this limit into 
the problem only in the combination 9/e 1 ^ 2 , therefore 
also the critical misorientation scales as 9 C ~ e 1 / 2 . 

As shown in Fig. |8j the range of repulsion is very 
narrow in PFC for e = 0.0923, which corresponds to 
e 1/2 = 0.304, differing from MD as seen in Fig. |] This 
result can be interpreted through the ratio of jsi/'fgb- 
A system with a larger ratio prefers grain boundaries to 
solid-liquid interfaces, narrowing the premelting range, 
while smaller ratios increase the range of misorientations 
that premelt. With this understanding, in order to obtain 
the same repulsive range of misorientations as exhibited 
in MD, we have to adjust e to decrease the PFC ratio to 
that of MD. Based on the analysis above, the material 
parameters for MD and PFC for e = 0.0923 are sum- 
marized in Table |l] These values are calculated using the 
procedure put forth in Ref.|3Sl exhibiting the discrepancy 
between PFC and MD. In PFC -y sl /B = 2.67 whereas in 
MD jsi/B — 1.59, where B denotes the bulk modulus. 



Based on the analysis in Eq. ( 62 ) above we can adjust this 
ratio in PFC by choosing e = 0.0329. In the dimensional 
units this is achieved by renormalizing S(qo) = 5.04 and 
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Quantity 


PFC bcc 


AE bcc 


MD bcc 


Cn (GPa) 


90.0 


90.0 


128.0 


C12 (GPa) 


45.0 


45.0 


103.4 


C 44 (GPa) 


45.0 


45.0 


63.9 


Bulk Modulus (GPa) 


60.0 


60.0 


111.6 


(Cn + 2<7 M )/3 








7 si (mJ/m 2 ) 


160.47 


144.26 


177.0 



TABLE I. Comparison of elastic constants and solid-liquid 
interfacial energies at the melting point predicted by the 
P FC model for e = 0.0923, AE and MD simulations for bcc 
Fe j44 | 45 | 54 | 82 | Note that thfi elastic consta nts for PFC and AE 

are the same. 



Quantity 


Original 


Rescaled 


Unit 


n 


0.0765 


0.0765 


A" 3 


u s 


0.72 


0.72 




qo 


2.985 


2.985 


A" 1 


a 


-2.136 


-1.274 


eVA 3 


X 


0.291 


0.488 


eVA 7 


9 


9.705 


258.646 


eVA 9 


e 


0.0923 


0.0329 




S(qo) 


3.01 


5.04 




C"(q ) 


-10.40 


-17.43 


A 2 


L 


1.968 • 10 9 


1.968 -10 9 


J/m 3 



TABLE II. Calculated constants and MD input parameters 
from the MH(SA) 2 potential. The Original values corre- 
spond to e = 0.0923 while the Rescaled values correspond to 
e = 0.0329. Calculated constants are derived in Appendix [C] 



C"(qa) = —17.43 A in comparison to the previous val- 
ues 3.01 and —10.4 A respectively, thereby leaving their 
ratio constant, which keeps the value of 7^; unchanged, 
see Eq. ((60 1. 



As seen in Fig. [10} for the new value of e the critical an- 
gle has shifted to a value of 8 ss 15°. The corresponding 
results for AE are shown in Figs. (TT( and [T2] Since the 
value of e is smaller here, the agreement between AE and 
PFC is significantly better than for the higher e value. A 
direct comparison is shown in Fig. |13[ The logarithmic 
plot Fig. [l4] shows that the decay ranges are very much 
comparable for AE and PFC. As demonstrated in Ref.1721 
the decay length also agrees with analytical predictions. 
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FIG. 9. (Color online) Liquid layer width as a function of 
grain boundary width for various misorientations and e = 
0.0329, obtained from PFC simulations. There is a transition 
from attractive to repulsive disjoining potentials at approxi- 
mately 6 = 15°. 
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FIG. 10. (Color online) Disjoining potential as a function 
of grain boundary width for various misorientations and e = 
0.0329, obtained from PFC simulations. There is a transition 
from attractive to repulsive disjoining potentials at approxi- 
mately e = 15°. 
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In Fig. [15] we also present three different MD disjoining 
potentials, over a similar range of angles as PFC and AE. 
We see that the range of of the interaction for PFC and 
AE is now comparable to the range of interaction in MD. 



FIG. 11. (Color online) Melt layer thickness for e = 0.0329, 
as computed from AE. There is a transition from attractive 
to repulsive disjoining potentials at approximately 8 = 14°. 
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FIG. 12. (Color online) Disjoining potential for e = 0.0329, 
as computed from AE. There is a transition from attractive 
to repulsive disjoining potentials at approximately 6 = 14°. 



FIG. 15. (Color online) Three different MD disjoining po- 
tentials. The points are the simulation data, the lines are 
exponential fits. The results are compared to PFC data for 
the rescaled value of e = 0.0329. 
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FIG. 13. (Color online) Disjoining potential comparison be- 
tween PFC and AE for e = 0.0329 for selected misorientations. 
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FIG. 16. (Color online) Histograms of PFC with noise grain- 
boundary liquid layer width as a function of three different 
undercoolings for a misorientation of 43.6°, e = 0.0923, and a 
cutoff of the noise at the scale of one unit cell (X m in ~ a). 
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FIG. 14. (Color online) Disjoining potential comparison for 
e = 0.0329 for selected misorientations. The exponential de- 
cay for the corresponding angles is similar for PFC and AE. 



C. Fluctuation effects 

As seen in Fig. [7] there is a profound difference in the 
range of the disjoining potential for MD and PFC with 
e = 0.0923. Above, we attempted to fix this difference 
through a rescaling of the parameter e. Another possi- 
ble route is to look at the role of fluctuations in PFC. 
Intuitively spoken, the underlying idea is that the ther- 
mal fluctuations lead to a broadening of the melt layer 
due to the asymmetric nature of the disjoining potential 
^(W). We quantitatively incorporate fluctuations into 
our calculation of the disjoining potential for one high 
angle boundary. As discussed above in Section |II B| we 
choose to cutoff the noise on the scale of one unit cell. In 
Fig. [16] histograms of the width distributions for three 
different undercoolings for a misorientation of 8 = 43.6° 
are shown. Using the procedure presented above in Sec- 
tion |IVB| we can extract the disjoining potential. Fig. 
17 shows a comparison of MD disjoining potentials, PFC 
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FIG. 17. (Color online) A comparison of PFC disjoining po- 
tentials with and without noise with an MD disjoining po- 
tential. The PFC results are for 9 = 43.6° and e = 0.0923 
for two different values of X m i n , while the MD results are for 
9 = 46.4° and are for three different undercoolings, a is the 
length of one lattice spacing. 



without noise, and noise based PFC for two different val- 
ues of Xmin- Akin to renormalizing e, the fluctuation 
effects clearly increase the range of the disjoining poten- 
tial. However, the strength of this renormalization is 
dependent on \ m %n- 

The introduction of fluctuations has the additional ef- 
fect of renormalizing the equilibrium chemical potential 
(// e g). Without noise e = 0.0923 and fi eq w -0.19609. 
The renormalized value of the equilibrium chemical po- 
tential is calculated by running cyrstal-liquid coexistence 
simulations, and we determine that for X m i n = la, 
Me<? ~ —0.18721 and for X m i„ = 2a, [i eq ps —0.19565, 
where a is one lattice spacing. 



D. Summary 



All the essential data are summarized in Figs. 18 and 
[T9| The plots show that the rescaling of e has a similar 
effect as the inclusion of thermal fluctuations, bringing 
the PFC and the MD results closer together. The log- 
arithmic plot shows that the interaction is indeed short 
ranged and decays exponentially for all simulation tech- 
niques. Nevertheless, there is still a significant difference 
in the slopes, i.e. the range of the interaction, between 
MD and the deterministic calculations. In comparison, 
the simulations with fluctuations exhibit a very similar 
behavior for short distances, because there the thermal 
fluctuations cannot compete with the overall strength of 
the repulsion. At larger grain separations, W > 1.5 nm, 
the data scatter becomes larger, and the change of slope 
may also be due to insufficient knowledge of the solid- 
liquid interfacial energy. Nevertheless, there is a clear 
tendency that (i) by better matching of the elastic con- 
stants, and (ii) incorporation of thermal fluctuations a 
better agreement of continuum methods and MD simu- 



FIG. 18. (Color online) Comparison of results for PFC and 
AE for the two different values of e with the PFC results with 
noise and the MD data. Additionally, the AE data for two 
crystals, which do not have a misorientation but which are 
just shifted against each other by half a lattice unit is shown. 
This curve is for a normal direction 21.8° off [100], which is 
the same as for the symmetric tilt grain boundaries, but here 
both grains are rotated in the same direction. 
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FIG. 19. (Color online) The same as Fig. 18 but on a log 



arithmic scale, to show the interaction range for the different 
models. 



lation results can be achieved. 

Finally, the plot also demonstrates that at on the rel- 
evant lengthscales the structural short range forces are 
substantially larger than London forces. Non-retarded 
dispersion forces between planar interfaces are asymp- 
totically described by the interaction potential 



V d 



H 



12ttW 2 ' 



(68) 



with the Hamaker constant H. For a characteristic value 
H = 10~ 21 J (see Ref.183]) the dispersion forces are of the 
order Vd ~ 10~ 2 mJ/m 2 at distances W ~ 1 nm, showing 
that their contribution is negligible. We also note that 
additional repulsive entropic forces are completely negli- 
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gible in comparison to the short-range structural forces 
computed here, see also the discussion in Ref. 1211 



VI. GRAIN BOUNDARY SHEARING 

In this section we use PFC to explore the relationship 
between premelting and the response of grain boundaries 
to shear stress. Dry symmetric tilt grain boundaries, or 
in PFC boundaries which are sufficiently undercooled, 
are known to follow a strict geometrical relationship be- 
tween the velocity of the applied shear parallel t o the 
boundary (tin) and its motion normal to itself (tvj 32 * 33 - 
This geometric model of coupling predicts two different 
branches of motion, /3^ 100 ) and /3(no) where 



V\\ = Pv n , 

£(100) = 2tan ( 
/3(iio) 



2 tan 

4 2 



(69) 
(70) 

(71) 



where 9 is the angle of total misorientation. Which cou- 
pling branch is preferred depends on the Burgers vector 
content of the grain boundary. For bcc low angle mis- 
orientaions with normal near (100), b = [100]a where 
a is the lattice spacing. These boundaries couple with 
j3 — /3(ioo)- F° r l° w angle misorientations with normal 
near (110), b = [110]a and (3 = /3(no)- For misorien- 
tations where the dislocations overlap, either coupling 
branch is accessible. In order to confirm the coupling 
branches for PFC we shear a bcc crystal with a symmet- 
ric tilt grain boundary using the dynamics presented in 
Eq. (13 1. These dynamics quickly relax elastic strain al- 



lowing us to shear at a faster velocity without deforming 
the crystal. To shear the crystal we add a term to the 
free-energy functional of the form, 

F' = F + J dVg(x) (ij-ij ) 2 = J dVf, (72) 

where g{x) is a normalized Gaussian of the form 



9{x) 



1 

V27TO- 2 



with a = a (one unit cell) and where 



(73) 



The expression %pf is calculated using Eq. (26 1. The 
added term to the free-energy effectively drags a strip 
of the crystal with a constant velocity of ±v. In our sys- 
tems we place one Gaussian strip to the right of the grain 
boundary shearing downwards, and one strip to the left of 
the grain boundary shearing upwards. We perform sim- 
ulations for a full range of misorientations 9 with renor- 
malized e — 0.0329 in order to better compare to MD re- 
sults and at an undercooling which ensures the physical 
contact of the two sides of the grain boundary (meaning 




FIG. 20. (Color online) An example of /J(ioo) symmetric 
tilt grain boundary coupling. The grain on the right is being 
sheared downwards with velocity v — —0.00025, while the 
grain on the left is being sheared upwards with velocity v = 
0.00025. The green atom starts in the right crystal, but due 
to the rightward motion of the boundary becomes attached 
to both crystals. It will eventually follow the pink atom and 
join a (100) plane of the left crystal. 
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FIG. 21. (Color online) Observed /3 as a function of misorien- 
tation 9 for e = 0.0329, and the ideal /3 factors for (100) 
and (110) coupling. These simulations were conducted at 
T/Tm ~ 0.7, and a net shearing velocity of U| = 0.0005. 



a negligible liquid layer 
vii = 0.0005 



and a net shearing velocity of 



In Fig. 



21 



we see that for low tempera- 
ture, PFC does reproduce the predicted ideal coupling, 
switching between the two branches at 9 = 45°. 

A recent studjES explored the existence of an alterna- 
tive grain boundary structure for low angle misorienta- 
tions with normal near (110). This "unpaired" structure 
splits the [022] a dislocations into pairs of [111] a/2 and 
[lll]a/2 dislocations. Simulations performed with this 
unpaired grain boundary structure exhibited the same 
observed coupling factor as the paired grain boundary. 
This is easily understood as the net Burgers vector con- 
tent of the two boundaries is the same. It is also possible 
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FIG. 22. (Color online) Observed /3 as a function of misori- 
entation 8 and homologous temperature for e = 0.0329. We 
shear with a net velocity Vu = 0.0005. 




FIG. 23. (Color online) The ratio of j3id ea i to the observed 
j3 as a function of misorientation for an undercooling of 0.46 
K corresponding to the top row of Fig. |22| Simulations are 
performed with vu = 0.0005 and e — 0.0329. There is a 
transition from ideal coupling to something approaching pure 
sliding as a function of misorientation. 



at temperatures closer to Tm, where the atoms on either 
side of the grain boundary are no longer in immediate 
contact, to observe behaviors other then pure coupling. 
An MD study of Ni (Ref. observed pure sliding for 
some misorientations as the melting temperature was ap- 
proached from below. In order to delineate these regions 
we perform a series of simulations over the full range of 
misorientations, and over a temperature range extending 
from an undercooling where we observe pure coupling at 
all 6 up to temperatures 1 K below Tm- 

We see in Fig. [22] that for temperatures closer to the 
melting temperature we observe coupling modes other 
than ideal coupling. There is a "v" shaped region where 
the observed j3 is larger than Pideai- 



In Fig. 23 we examine the observed coupling factor 
close to Tm- For low angle misorientations, where the 
liquid layer width is smaller, we still observe pure cou- 
pling. High angle misorientations however, have coupling 



factors much larger than ideal. 



VII. CONCLUSIONS 

We have used the PFC and AE methods to investigate 
the equilibrium prcmclting and noncquilibrium shearing 
behaviors of [001] symmetric tilt grain boundaries (GBs) 
at high homologous temperature in classical models of 
bec Fe, and also compared our equilibrium results to MD 
simulations. 

At a qualitative level, our findings a re con sistent with 
the findings of previous PFCPand MEEUm studies. We 
find that the disjoining potential can be either repul- 
sive over an intermediate range of misorientation (f5 m i n < 
< $max), or exhibit an attractive minimum outside this 
range (9 < # min or max < < 90° ) . At a more quantita- 
tive level, we find that the repulsive range of misorienta- 
tion is not well predicted by the PFC model for the value 
of e matched to liquid structure factor properties. The 
range predicted by PFC is too narrow compared to MD^]. 
This disagreement stems from the fact that the PFC 
model is too simplified to reproduce, for the same value of 
e, the correct liquid structure factor properties, the cor- 
rect solid-liquid interfacial free-energy, and the correct 
elastic moduli. However, if e is chosen to have a lower 
value so as to match the ratio of bulk modulus and j s i in 
the MD simulations, the PFC model predicts accurately 
the range of misorientation where GB exhibit a diverging 
layer width at the melting point. Since the elastic mod- 
uli are a major determinant of the GB free-energy at the 
melting point, this result indicates that elastic properties 
and the solid-liquid interfacial free-energy play a domi- 
nant role in setting the range of misorientation where 
GBs premelt. In the AE framework, the range of misori- 
entation where GBs remain dry, with a minimum in the 
disjoining potential, is predicted to shrink to zero in the 
small e limit with (9 m j n and 90° — max both scaling as 

Furthermore, we have found that with the lower value 
of e, which predicts the correct premelted range of mis- 
orientation, the PFC model reproduces the expected GB 
shearing behavior as a function of homologous tempera- 
ture and misorientatiorPS. Namely, for low enough ho- 
mologous temperature, GBs exhibit pure coupled mo- 
tion to a shear stress with a discontinuous change of the 
coupling factor as a function of 9 that reflects a tran- 
sition between two coupling modes. In contrast, for 
high enough homologous temperature (T/Tm > 0.8 in 
the PFC model), partial sliding occurs, thereby causing 
the coupling factor f3 to be lower than its ideal value. 
The range of misorientation where /3 < fiideai increases 
with temperature but remains finite even at the melting 
point since low angle GBs with well separated disloca- 
tions exhibit pure coupled motion even when the dislo- 
cation cores are partially melted. 

Those results support the view that, for symmetric tilt 
GBs, partial sliding is primarily due to the formation of a 
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premelted layer at high homologous temperature. How- 
ever a recent combined PFC and MD study^ has shown 
that asymmetrical tilt GBs can exhibit partial sliding 
even for low temperatures where such a layer is absent. 
Therefore further studies remain needed to obtain a more 
complete understanding of GB sliding mechanisms as a 
function of temperature and GB bi-crystallography. 
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Appendix A: Integration schemes 
1. PFC with noise 

The evolution equations for non-conserved dynamics 
with noise are given in Eq. ^9b. Due to the high order 
derivative terms, it is computationally more tractable to 
conduct the simulations in reciprocal space. We apply 
a Fourier transform as in the deterministic dynamics in 
the Appendix of Ref. [29l 



d t ip k = L k ipk + 9k> 
L k = (e-l) + 2fc 2 -fc 4 , 

= fk + Vk: 



(Al) 
(A2) 



(A3) 
(A4) 



We can now proceed as in Ref. [29] 

y] fe = u{t)e L "\ (A5) 
d t i) k = L k e Lkt u{t)+g k , (A6) 
d t u(t) = e- L ^~g k . (A7) 

Integrating u{t) in time gives us the expression, 
*t+At 

u(t + At) - u(t) = J d£ exp(-L k t)g k (t'), (A8) 

where f k (t') can be expanded around t! — t. Since the 
term fj k is an instantaneous value over the interval At 
it can be treated as a constant during integration. Our 



final expression for the dynamics in a continuous case is 
then 



Mt + At) 



,L k (t+At) 



t+At 



dt't 



-LiA 



/*(*) 



f k (t)~fk(t-At) 
At 



e L k At^ t)+ f_]M^L k At 

L k 

f k (t)-fk(t-At) 



(t' -t) + fj k 
1) 



AtL\ 



L k At 



- 1 - AtLi 



1 



L k At 



-Vk- 



(A9) 



Special care must be taken to properly normalize the 
noise in spectral space. In real space the average noise 
and two point correlation function of the noise are 



(A10) 



2 J^5(f-r*W-t'). (All) 



The last equation includes the strength of the Gaussian 
noise (2fcs? 1 ) converted into PFC dimensionless units. 
While these equations accurately describe the addition 
of noise for a continuous system in real space, care must 
be taken to normalize the noise properly in a discretized 
spectral space, specifically for the fast Fourier transform 
which we use in this study. We start for simplicity in one 
dimension where 

(r,(x, t)r,'{x>, 0) - 2k ^5(x - x')S(t - t'). (A12) 

For a discretized system we write 

2k B Tg Sji S nm 



\ 2 q r ^ Ax At 



(A13) 



where j,l,m,n are integers, x{x') = j(l)Ax, and t(t') = 
n(m)At. We then define the discretized Fourier trans- 
form of the noise to be 



N-l 

E 

N-l 



e -2irijr/N^n 



1=0 



(A14) 



where N is the number of points in the system and r(s) 
is an integer that takes on all values from —N/2 + 1 to 
N/2 inclusive, r can be related to the frequency through 
the formula k = 2irr /N Ax. The noise in Fourier space 
must be normalized separately for both real and imagi- 
nary parts. 



iri"" 1 



(A15) 
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where 



Vr"" = E C " S 
3=0 
N-l 

^ - E 

3=0 

m?) = 



N r > » 



■ 2irjr „ 



(A16) 
(A17) 

(A18) 



It can easily be seen that the two crossterms do not con- 
tribute to the overall noise, 



(v. 



2k B Tg 5 nm ^-^ Sji 2njr . 2nls 
x cos sin - 



A 2 «5 At ^ Ax ~~~ N 
y ° 3,1=0 



N 

2'KJS 



2k B Tg S nm ^-^ 2njr . 

: ^— c : : > COS Sin 

\ 2 qk AtAx ^ N N 
H ° j=o 



(A19) 



(A20) 



which is equal to zero for all integers r and s. The other 
two terms do contribute to the noise, 



2k B Tg 5, 



E 



2-njr 2irls 
t s / t— cos „ cos ■ 

A 2 g£ At ^ Az 



iV 



x cos — cos 



3=0 



\ 2 q% AtAx 

2k B Tg5 mn ^ 
X 2 q^AxAt' 



N 



N 
2njs 



N N , 

~1? + -W{°rQ + <V JV/2. 



(A21) 

(A22) 
(A23) 



and 



~n,I ~m,I\ 



2k B Tg S mn ^-^ Sji . 2njr . 2-rtls 

Z ^ n r; 7 > ~T — Sill — Sill 

\ 2 ql At 4< Ax 

2k B Tg S nm 
\ 2 ql AtAx 



3,1=0 
N-l 



N 



N 



E 



. 27rj'r . 2njs 
sin — — — sm ■ 



3=0 



N 



N 



2k B TgS 

II 

X 2 q^AxAt' 



N N 



Srs I ~Z ^(^rO + fir N/2) 



(A24) 

(A25) 
(A26) 



As can be seen both the imaginary and real terms con- 
tribute equally for r, s 7^ 0, but the r, s = and 
r, s = N/2 terms have the full magnitude only on the 
real part of the dynamics, and are zero for the imagi- 
nary part. The logic extends easily to three dimensions 
with N/Ax replaced by (N x N y N z )/(AxAyAz). As pre- 
sented in section |II B| in order to ensure that the noise 
based contribution to the free-energy is less than the free- 
energy of a noiseless system we must cutoff the noise at 
a length scale close to one lattice spacing. We choose the 
following form for our cutoff 



tip* 




{5r0 + firN/2)^ X 



(A27) 




(S r0 + 5 rN / 2 ) j x 



(A28) 

where £ controls the width of the cutoff, A mi „ is the min- 
imum noise wavelength, T = 2k B Tg / (\ 2 q$), and G are 
Gaussian random numbers with a standard deviation of 



2. Wave Dynamics 

Here we present the algorithm used to implement the 
modified PFC dynamics introduced in Ref. [721 We re- 
peat the free-energy and dynamical equation for locally 
conserved dynamics 



F 



/ 



dr 



i { - e + { v 2 + i) 2 )ip + ^ 



v> 2 

~2 



d tt i> + ad t ip = V : 



,6F 
6i/>' 



(A29) 
(A30) 



As in Ref. [2H1 we avoid calculating the gradient terms in 
real space by solving the equation of motion in Fourier 
space. We define a Fourier transform 



(A31) 



We multiply both sides of the dynamical equation by 
exp(ifc • r) and integrate. The equation of motion in 
Fourier space is now 



duipk + otd t ipk = Lk^k + fk, 



where 



fk 



fc 6 + 2fc 4 + (e-l)fc 2 , 



(A32) 

(A33) 
(A34) 



Following the methodology of Section 3 of Chapter II of 
Ref. I84| we now choose to rewrite the dynamics as 



V^fc = u k , 

u k = -au k + L k ip k + f k . 



(A35) 
(A36) 
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We now solve the homogenous form of Eq. ( A32 ) by writ- 
ing 



where 



Mi = -a/2 + ^a 2 /4 + L k , 
fi 2 = -a/2-Ja 2 /A + L k , 



(A38) 
(A39) 



and where cti and a 2 are functions of time and restricted 
to satisfy 



dai 
~dt' 



.An* 



da 2 
~dt 



3"»* - 0. 



(A40) 



Combining Eq. ( A37 1 and Eq. ( A32 ) we derive the rela- 
tion 



da i 



(A41) 



Additionally combining Eq. (A41) and Eq. (A40) and 



integrating over one time step we obtain 



1 



Al 



«.,. = - / e-^ f k (t')dt' + a w , (A42) 

Mi - M2 Jo 

At 

e-^ t 'f k (t , )dt' + a 2 o, (A43) 



-1 



« 2 



Mi - M2 Jo 



where etio and 020 are integration constants. Now com- 
bining everything 



i) k = ai {t)e^ 1 + a 2 (t)e^ t , (A37) $ k {t + At) 



1 



Mi - M2 



/t+At 
e" Mlt ' 'h{t')dt' 



u k (t + At) 



/t+At 
e-» 2t ' f k {t')dt' 

+a w e^ t+At ^ + a 20 e^ t+At \ 
1 



(A44) 



Mi - M2 



/t+At 
e~ Mlt ' 'fk{t')dt' 



-M2e 



t+At 

/i 2 (t+At) / „-M2t 



e-^ 1 f k (t')dt' 



+^a 10 e^ t+A ^ + n 2 a 20 e^ t+At K 



(A45) 



We integrate over the range t = to t = At and require 
that at time t = 

4>k(t = 0) = i>ko = aw + a 2Q , (A46) 
u k (t = 0) = u fc0 = Mi fl io + M2O20- (A47) 

Solving for the integration constants aio and a 2 o 
(u k0 - iJL 2 i>ko) 



a 2 o 



Mi - M2 
-(uko ~ Mjjfeo) 
Mi - M2 



(A48) 
(A49) 



To calculate the integral terms we use the approach pre- 
sented in the Appendix of Ref.[Sn|by expanding f k around 
t'=t 



^(t+At) f t+At 



Mi - M2 Jt 



f k {t')dt' - 



^(t+At) f t+At 



Mi - M2 Jt 



dt'e-^' x f k (t) + 



A(t) -/*(*- At) 



At 



(t'-t) 



(A50) 



Integrating the first term in the integral 



pi(t+At) ft+At 



Mi - M2 Jt 



Mi(Mi - M2)* 



dt'e-^ f k (t) = ^ he^ t+At) = ^ ^ /*(*)• 

Mi (Mi - M2) 



(A51) 



Integrating the second term in the integral we see 



e , l( t+At) f t+At ,^(t)-A(t-At) 



Mi - M2 Jt 



df 



, M , ,„.,. _ = (f k (t) - f k (t - At)) f e^ At 1 1 

Mi(Mi-M2) \MiAt ^lAt 



Af 



(A52) 
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The full dynamics then take on the form 

4> k (t + At) = — 1 — (-(Mt)v2 - u k (tyK lAt + (Mt)m - u k (t))e^ At ) 
fit - fl 2 \ ' 



hit) f<JllM „ , (f k (t)-f k (t-At)) fe^ At 1 



Mi(Mi - Ma) Mi(Mi _ Ma) \ At /iiAt 

/*(*) , oMa At ^ (Mt)-f k {t-At))fe^ At 1 



(e ^ t _ 1} _ w»v " _ _ i 1 (A53) 



Ma(Mi-Ma) ^2(^1-^2) \M2Ai ^ 2 At 

u fe (< + At) = (-Hi{Mt)^2 - u k (t))^ lAt + l*»$k{t)Hi - Uk(t))e" 2At ) 

Ml — M2 \ ' 

+ AW (e , lA1 _ 1} + (A(t)-A(t-At)) f e^ At 1 A 
(M1-M2) (M1-M2) \MiAi Mi At / 

AW ,„ a At ^ (/*(*)-/*(*- At)) /e^ Ai 1 



(e ^ 2 At _ 1} _ ^ ^ (A54) 

(Mi _ Ma) (M1-M2) \fJ-2At (I? At 



The non-conserved dynamics take on a similar form 

(5F 



5^ 



/i, 



(A55) 



where is an externally imposed chemical potential. The 
only changes to the dynamics are in the operator L/. and 
the non linear term A 



L fc = -fc 4 + 2fc 2 + (e-l), 



(A56) 
(A57) 
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Appendix B: Liquid layer width for PFC with noise 

The inclusion of noise in the PFC model demands an 
alternative definition of the melt layer thickness W. The 
method of calculating the liquid layer width for PFC 
without noise by using the excess mass no longer works 
once we have added noise. This is due to the relatively 
large instantaneous fluctuations in the amount of mass 
in the system. Instead we introduce a wavelet transform 
with methodology similar to a transform used to measure 
the local orientation of a two-dimensional PFC hexagonal 
crystal 85 . In our case we extract an orientation indepen- 
dent field which differentiates continuously between solid 
and liquid. We define a convolution of the square of the 
gradient of the density field ip with a normalized Gaus- 
sian function, 



FIG. 24. Surface plot of the order parameter p(f) defined 



where 



p(r) = / G(\f-r*\)\Vip(T*)\ 2 (lr' 



G = 



(Bl) 



(B2) 



and (T2 — 2.5a where a is one lattice spacing. The con- 
volution returns an orientation independent field whose 



in Eq. ( B2 \ through a solid-liquid interface for PFC without 
noise. The graph shows the order parameter plotted for a 
single slice in the x-y plane (for a single value of z). In the 
solid we see p s ~ 1.22 while in the liquid pi — 0, with a tanh 
like interpolation between the two values, a is one lattice 
spacing. 



amplitude depends on how close to a perfect solid the 
crystal is. As seen in Fig. [24] this results in a smooth 
transition from the solid to the liquid through a solid- 
liquid interface in the deterministic case. The behavior 
of the order parameter is different for a grain-boundary 
interface in the deterministic case, as seen in Fig. [25] p 
never reaches its minimum value of in between the two 
solid regions. This behavior is consistent with an inter- 
face where the two grains are not well separated. 

The behavior of the order parameter when including 
noise is more complicated, as seen in Fig. [26] Instead 
of attempting to determine the liquid layer width W in 
one calculation, we separate the order parameter into an 
array of functions p(x,yo, zq), where x takes on its full 
range of values and y and Zq refer to a particular point 
in the y — z plane. We then fit the following equation 
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FIG. 25. Surface plot of the order parameter p(f) defined in 
Eq. ( |B2[ ) through a grain-boundary interface for PFC without 
noise. The graph shows the order parameter plotted for a 
single slice in the x-y plane (for a single value of z). In the 
solid we see p~ s « 1.22, however the order parameter does not 
take on its minimum value of p = between the two crystals. 
This is consistent with an interface where the two grains are 
not well separated, a is one lattice spacing. 
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FIG. 26. Surface plot of the order parameter p(r) defined 
in Eq. ( |B2[ ) through a grain-boundary interface for PFC with 
noise. The graph shows the order parameter plotted for a sin- 
gle slice in the x-y plane (for a single value of z) . As compared 
to the deterministic case, there is more variation in the order 
parameter in the solid and the central liquid like region. This 
makes it necessary to average over several points to identify 
the liquid and solid values of the order parameter, as detailed 
in the text, a is one lattice spacing. 



pi(yo,z ) 



Ps(yo,z ) + pi{y 0l z ) 
2 

Ps(yo,z ) - pi(y ,z ) 



(1 - ta,nh(bi(x - x t ))) 



(1 + tanh(6 r (x — x r ))), 
(B3) 

to p(x, yo,zo) for each combination of yo and zg. bi, b r , xi, 
and fit parameters while p s {yo, zq) is the value of 

the order parameter in the solid and is obtained by aver- 
aging p{x, yo; Zq) over 5 different values of x, where x is 



several lattice spacings away from the interface. pi{yo, zq) 
is defined as the average of the two smallest values of 
p(x,y ,z ). 

W(yo, zq), the liquid layer width for each set yo, zq, is 
defined 9iS X f X l • The total liquid layer width W is then 
the average of W(yo, zq) over all yo, zq. 



Appendix C: Linking PFC and DFT 

To obtain the value of the only adjustable parameter, e, 
we derive an expression to relate e to material parameters 
taken from MD simulations, namely the peak and curva- 
ture at the peak of the liquid structure factor. While 
already derived in Ref . [??] we use method presented in 
Ref. H5] for the two-mode PFC model tha t mor e directly 
relates these quantities to classical DFTS^3186| We var y 
ip around its liquid value ^ = ipi + Sip , where ipi is the 
density in the liquid, as determined from the phase co- 
existence conditions between solid and liquid, and sub- 
stitute into Eq. (8j). We see that (dropping terms of Sip 
higher than quadratic order) 



AT - 



9 



dr 



Sip 



[a + 3^ 2 A % 4 + A(V 2 + ql) 2 ]SiP 



Defining the Fourier transform, 

sip • d% 



Sipke 



ik-f 



(2tt) 3 / 2 

and substituting this definition into Eq. (CI I 
A<7q f f dkdk' SipkSipk' ' 



(CI) 
(C2) 



A J 7 



9 J J {2nf 2 
+g 2 ) 2 ) / dfe^ + ^ p 



(a + 3tpfXq^ + X(-k A 



Ag 4 



dk 



Sip k Sip- 



<?(2tt) 3 / 2 
+A(-fc 2 + 9o 2 ) 2 ] 



[a + 3iPf\q^ 



(C3) 



We can compare this equation to a standard formulation 
from classical DFT, 



DFT 



k B T 
2 

Sn(r) 



dfdr 

Sir-r 1 ) 
n 



C(\f-7\) 



Sn{f), 



(C4) 

where the density variations 8n(r) are related to the PFC 
order parameter via Eq. (|6| 



Sn(r) = n(r) — n = S(p(f) 



'Ago 4 



Sip(r 



and 



C(k) = n / dfC(\f\)e 



-ik-; 



(C5) 



(C6) 
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is the Fourier transform of the direct correlation function. 
Substituting in the Fourier transforms as before 



DFT — 



g 2n (2vr)3/2 



dk5i> k 5^_ k [l-C{k)]. 



Equating the two free-energy densities and using the ex- 



pression for the liquid structure factor S(k) 
C(k)] we see that 



S{k) 



k R T 



n {a + 3\qltf + \(-k 2 + q 2 ) 2 }- 



Evaluated at the peak of the structure factor (k 
see, 



a + 3 Agg 



k B T 
noS(qo) 



1/[1- 

(C8) 
■ g ) we 

(C9) 



The multiscale expansion using e 1 / 2 as small parameter 
is discussed ir!2Q in detail, and we repeat only the results 
here. To this end the density field is tp in expanded in 
powers of e 1 / 2 , 



and we can expand the average solid and liquid densities 
in the same way 



$1 



^ s0 e 1/2 + ^ sl e + ^ s2 e 3/2 + 
^ w e 1/2 + ^i 1 e + ^i 2 e 3/2 + 



(Cll) 
(C12) 



From the common tangent construction one obtains at 
the different powers of e 



Ipsl = Ipll = 0. 



(C13) 
(C14) 



We see that in the small e limit the size of the solid- 
liquid coexistence region A-ip = ^ s — ipi « (ip S 2 — ipi2)e 3 ^ 2 
is much smaller than the average density which scales like 



e 1 / 2 justifying our dropping of higher order terms in the 
small e regime. 

Using Eq. Q to substitute for a and the small e 
justification for dropping higher order terms to define 
$1 = ip c e 1/2 we get, 



-k B T 



noS{q )\qo(l - 3i/> 2 



(C15) 



Taking the second derivative of C(k) and evaluating at 
k = qo we obtain a value for A, 



A 



k B TC"(q ) 



2o n o 



(C16) 



see also Eq. (40 1, and we thus get 



6 (l-3^ 2 )q 2 S(qo)C"( qo )\ JP 17) 
which matches an expression originally derived irP^. Us- 
ing the values for ip c , 5 (go) an d C"(qo) reported irP31(see 
also Table |TTJ) we get e = 0.0923. 

In order to convert the dimensionless PFC results back 
into dimensional units we must also define a, A and g. 



(CIO) a and A can be defined through Eqs. (C9| and (C16) 



respectively. To calculate g we first recognize that 



noUi 



/2 A 



(C18) 



where m are the dimensional amplitudes (as opposed to 
the dimensionless amplitudes Ai) and A° = e~ l / 2 A s . 
Substituting — u s , which can be obtained from MD 
simulations^ we obtain, 



.9 = Ag 4 £ (-^ c+ ±y/K=uj%j /(n 2 u 2 ). (C19) 

This is a correction from the expression first presented 
irP. 
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